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The structure of symplectic integrators up to fourth-order can be completely and analytical un- 
derstood when the factorization (split) coefhcents are related linearly but with a uniform nonlinear 
proportional factor. The analytic form of these extended-linear symplectic integrators greatly simpli- 
fied proofs of their general properties and allowed easy construction of both forward and non-forward 
fourth-order algorithms with arbitrary number of operators. Most fourth-order forward integrators 
can now be derived analytically from this extended-linear formulation without the use of symbolic 
algebra. 



I. INTRODUCTION 

Evolution equations of the form 



w(t + e) =c £{T+v) w(t), (1.1) 



=1 



via a well chosen set of factorization (or split) coefficients {ti} and {vt}. The resulting algorithm is then nth order 
because the algorithm's Hamiltonian is T + V + 0(e n ). By understanding this single approximation, computational 
problems in diverse fields of physics can all be solved by applying the same algorithm. 

Classically, (|1.2f> results in a class of composed or factorized symplectic integrators. While the conditions on {ti} and 
{vi} for producing an nth order algorithm can be stated, these order conditions are highly nonlinear and analytically 



where T and V are non-commuting operators, are fundamental to all fields of physics ranging from classical 
mechanic9i*2i2i^, electrodynamics^, statistical mechanics&SiiSiii to quantum mechanicsiSii&ii. All can be solved 
O^l ■ by approximating e £ ^ T+v ) to the (n + l)th order in the product form 
> 

g ! e z(T+V) = Y[ e U*T e v<eV + ( £ n+lj (L2) 

in : 
o ■ 

■ opaque. In many casesi^iiSiiSiii, elaborate symbolic mathematical programs are needed to produce even fairly low 
order algorithms if N is large. In this work, we show that the structure of most fourth-order algorithms, including 

S nearly all known forward ({ti, Vi} > 0) integrators, can be understood and derived on the basis that {vi} and {ti} are 
linearly related but with a uniform nonlinear proportional factor. This class of extended-linear integrator is sufficiently 
^ , complex to be respresentative of symplectic algorithms in general, but its transparent structure makes it invaluable 
for constructing integrators up to the fourth-order. In this work we prove three important theorems on the basis 
of which, many families of fourth-order algorithms can be derived with analytically known coefficients, including all 

■ known forward integrators up to N — 4. 

_ _ i 

II. THE ERROR COEFFICIENTS 

The product form l|1.2|) has the general expansion 

N , 

Y[ e tieT e VieV = exp I ee T T + ee v V + e 2 e TV [T, V] 
i=i ^ 

+ e 3 e TT v[T, [T, V}} + e 3 e VTV [V, [T, V]] + ■ (2.1) 

We have previously^ described in detail how the error coefficients e-r, ey, ctv , exTV, and eyrv can be determined 
from {U} and {vi}: 

N N 

e T = y^ j t i , ey^^Vi, (2.2) 

i=l i=l 
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(2.3) 



i=i 
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— + -e T v + e T rv 



1 1 

T^j + 2 eTV ~ evTV 



1 N 



i=l 
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(2.4) 
(2.5) 



where we have defined useful variables 



and the backward finite differences 



with property 



i N 

si = tj , Ui = Vj , 
i=i 3=t 



AT 

E 



(2.6) 



(2.7) 
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We will always assume that the primary constraint ey = 1 and ey = 1 are satisfied so that l|2.8|) sums to unity. 
Satisfying these two primary constriants is sufficient to produce a first order algorithm. For a second-order algorithm, 
one must additionally forces erv — 0. For a third-order algorithm, one further requires that errv — and evrv = 0. 
For a fourth-order algorithm, one has to satisfy the third-order constraints with coefficients ti that are left-right 
symmetric. (The symmetry for Vi will follow and need not be imposed a priori.) Once the primary conditions 
Bt = 1 and ey = 1 are imposed, the constraints equations (|2.3|I - H2.5|I are highly nonlinear and difficult to decipher 
analytically. In this work, we will show that l|2.3|l can be satisfied for all N by having {v{} linearly related to {ti} (or 
vice versa). The coefficients exTV an d evTV can then be evaluated simply in terms of {ti} (or {v^}) alone. This then 
completely determines the structure of third and fourth-order algorithms. 



III. THE EXTENDED-LINEAR FORMULATION 



The constraint erv = is satisfied if 



N 1 

E VsiUi = o- 



(3.1) 



If we view {ti} as given, this is a linear equation for {ui}. Knowing Q2.8|) . a general solution for Ui in terms of Sj and 
Si-i is 



M V7«™ M 1 



n=l 



(3.2) 



The coefficients C n respresent the intrinsic freedom in {vi} to satisfy any constraint as expressed through its rela- 
tionship to {ti}. The expansion l|3.2|) is in increasing powers of and Sj_i. If we truncated the expansion at M = 2, 
then for i ^ 1, m is linearly related to {s^}, i.e., 



Vs 2 

^ = Ci + = Ci + C 2 {si + Si-!). 

For i = 1, since we must satisfy the primary constraint ey = 1, we must have 

Ml = 1. 



(3.3) 



(3.4) 
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In this case, the constriant (|3.1|) takes the form 

N 1 

5^Vai«i = ti + Ci(l-ti) + C a (l-t?) = 2- (3-5) 

i=l 

The complication introduced by u\ = 1, in this, and in other similar sums, can be avoided without any loss of 
generality by decreeing 

tt = 0, (3.6) 

so that (|3.5(l remains 

C 1+ C 2 = 1 -. (3.7) 

For i ^ 1 ^ N, lO implies that 

Vi = -C2(ti+t i+1 ). (3.8) 
Since v\ = u\ — u 2 = 1 — C\ — C 2 t 2 , by virtue of (|3.7|l . 

«i = ^ + C , 2 (l-t 2 ). (3.9) 
Similarly, since vn = % = C\ + C 2 (2 — tjy), we also have 

ww = ^ + C 2 (l-tjv). (3.10) 

Given {ti} such that t\ = 0, the set of {vi} defined by l|3.8|l - H3.10|l automatically satisfies ey = 1 and e^y = 0. If C 2 
were a real constant, then {v t } is linearly related to {ti}. However, in most cases C 2 will be a function of {ti} and 
the actual dependence is nonlinear. But the nonlinearity is restricted to C 2 , which is the same for all Vi. We will call 
this special form of dependence of Vi on {ti}, extended-linear. For a given set of ti, (|3.8() - (|3.10() defines our class of 
extended- linear integrators with one remaining parameter C 2 . 

For extended-linear integrators as described above, one can easily check that the sums in (|2.4Jl and (|2.5|l can be 
evaluated as 

N 1 

Va?Uj = d + C 2 + gC 2 = -+ gC 2 , (3.11) 

i=l 

N 1 

V S . lMj 2 = (d + C 2 f + 9 Cf = - + gC 2 2 . (3.12) 

i=i 

Again the complication introduced by u\ = 1 is avoided by decreeing t\ = 0. The quantity g is a frequently occuring 
sum defined via 

i=l 



w i 

= SiSi-i(si - Si-i) = -(1 - (3-14) 



with explicit form 

N 

9 

i=i 

where 

AT 

fc = £*<- (3-15) 

i=l 

Much of the mechanics of dealing with these sums have been worked out in Ref^£. However, their use and interpre- 
tation here are very different. From (|2.4() and Q2.5|) . we have 

e T TV = ^ + ^.gC 2 , (3.16) 

evTV = ^-\gC 2 2 . (3.17) 

Both are now only functions of {ti} through g. 
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IV. FUNDAMENTAL THEOREMS 



We can now prove a number of important results: 
Theorem 1. For the class of extended-linear symplectic integrators defined by t\ — and (|3.8|) - (|3.1U|) . if {i,} > 
for i ^ 1 such that e-p = 1, then zttv 7^ &vtv- 
Proof: Setting exrv = ^vtv produces a quadratic equation for C 2 , 

Cl+C 2 + -^-=0 (4.1) 



12g 



whose discriminant 



D = b 2 - 4ac = — (4.2) 

1 - Sg 

is strictly negative (since if ey = 1, then 1 > Sg > 0). Hence no real solution exists for C 2 . This is a fundamental 
theorem about positive-coefficient factorizations. This was first proved generally in the context of symplectic corrector 
(or process) algorithms by by Chinii and by Blanes and CasasiS. If exTV can never equal evTVi then no second 
order algorithm with positive coefficients can be corrected beyond second order with the use of a corrector. 

As a corollary, for {Uyi} > 0, exrv and evrv cannot both vanish. This is the content of the Sheng-Suzuki 
TheoremSSiSi: there are no integrators of order greater than two of the form (|2.1H with only positive factorization 
coefficients. Our proof here is restricted to extended-linear integrators, but can be interpreted more generally as it is 
done in Refi^. Blanes and Casas^S, have also given a elementary proof of this using a very weak necessary condition. 
Here, for extended-linear integrators, we can be very precise in stating how both errv and evrv fail to vanish. We 
have, from l|3.16[) . if cttv = 0, then 

C2 = -2o 1 ^)' evTv= -YAj^6iy (43) 

Similarly, from l|3.17jl . if evrv — 0, then 

C ' 2 = -^T=r 9 > ZTTV = y 2 (l-V^). (4-4) 

Satisfying either condition forces C 2 to be a function of {ti} through Sg. From Refpi&, we have learned that the value 
given by (|4.3(l is actually an upperbound for evrv if {U>i} > and exrv — 0- Similarly, in general, the value given 
by H4.4fl is a lower bound for errv if {U>x} > and evrv = 0. Our class of extended-linear integrators are all 
algorithms that attain these bounds for positive £j>i. Note that in (|4.4() we have discarded the positive solution for 
C2 which would have led to negative values for the vt coefficients. 

For the study of forward integrators where one requires {ii>i} > 0, it is useful to state (|4.3(l as a theorem: 
Theorem 2a. For the class of extended-linear symplectic integrators defined by H3.8fl - 13.10fl with ti = 0, &t = 1, and 
C 2 , e VT v given by 



r - 1 

<^2 — — 777, evTV 



one has 

N 



f[ e tisT e^ v = exp (e{T + V)+ e 3 e VTV [V, [T, V}] + --^. 



(4.6) 



For t\ = 0, the first operator e VieV classically updates the velocity (momentum) variable. Theorem 2a completely 
described the structure of these velocity-type algorithms. 

If one now interchanges T «-> V and {ti} «-> {vi}, then [T, [T, V]} transforms into [V, [T,V]] with a sign change. 
Hence, one needs to interpret errv in 14-4(1 as — evrv, yielding: 
Theorem 2b. For the class of extended-linear symplectic integrators defined by 

h = - + C 2 {l-v 2 ), t N = - + C 2 (1 -v N ), ti = -C 2 (vi + Vi + i), (4.7) 
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with v± = 0, ey = 1, and 



^ = -^7, e V TV = -^(l-A (4.8) 



where 

JV 

one has 

]Je VieV e tiST = expfe{T + V)+ e 3 e VTV [V, [T, V}] + ■■ (4.10) 



For v-y = 0, the first operator e tlsT classically updates the position variable. Theorem 2b completely described the 
structure of these position-type algorithms. 

In both Theorem 2a and 2b, one obtains fourth-order forward algorithms by simply moving the commutator 
[V, [T, V]) term back to the left hand side and distribute it symmetrically among all the V operators^. 

If some ti were allowed to be negative, then both cttv and evrv can be zero for Sg = 0. For both (|4.3|l and 14.4(1 
we have 

C 2 = -~ (4.11) 

and 

Vi = ^(ti+t i+1 ). (4.12) 

The latter is now true even for i = 1 and i = N. This is not an coincident, from i|3.16[) and l|3.17[l . if we set C2 = —1/2, 
then 

e T TV = 2 e VT v = - 7 = J^ s 9- (4-13) 

Since C2 here is a true constant, {wi} is linearly related to {U}. We can formulate this explicitly as a theorem for 
negative-coefficient factorization yielding truly linear algorithms: 
Theorem 3: If 

Vi = ^(U + t i+l ) (4.14) 

such that ti = 0, then 

N / 3 \ 

II e^ T e^v = exp ee T (T + V)+ £ -Sg (2% [T, V}} + [V, [T, V]}) + ■ ■ ■ . (4.15) 
i=i ^ ' 

Both commutators now vanish simultaneously if <5<7 = 0. 

An immediate corollary is that if 5g were to vanish, then there must be at least one tk < such that t\ + tf. , x < 
or t\ +t 3 k _ 1 < 0. Since 

(a; 3 + y 3 ) = (x + y)[^y 2 + (x - ^y)\ 

x 3 + y 3 < x + y < 0. We therefore must have tk + ife+i < or tk + ife-i < 0. From l|4.14|l . this implies that Vk or 
Vk-i must be negative. Thus an algorithm of order greater than two of the form l|4.15|l must contain at least one pair 
of negative coefficients ti and vj. In its general context, this is the Goldman-Kaper result-. Our linear formulation 
here is more precise: if U is negative, then at least one of its adjacent Vi must be negative. If only one tk is negative, 
then both of its adjacent Vi must be negative. 
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V. THE STRUCTURE OF FORWARD INTEGRATORS 



Theorems 2a and 2b can be used to construct fourth-order forward algorithms with only positive factorization 
coefficients. These forward integrators are the only fourth-order factorized symplectic algorithms capable of integrating 
time-irreversible equation such as the Fokker-Planckifl^ or the imaginary time Schrodinger equation^2£i2&. Since 
it has been shown thati^ currently there are no practical ways of constructing sixth-order forward integrators, these 
fourth-order algorithms enjoy a unique status. 

For N = 3, for a fourth-order algorithm, we must require t 2 = t% = 1/2. Theorem 2a then implies that 

12 1 

Vi = v 3 = g, v 2 = -, and e VT v = —=^- (5.1) 

By moving the term s 3 evTv[V, [T,V]] back to the LHS of l|1.2fl and combined it with the central V, one recovers 
forward algorithm 4A 2 L2&. For N = 4 with t 2 = t 3 = t 4 — 1/3, we have 

13 1 , . 

vi = Vi = g, v 2 = v 3 = -, and e V rv = ^^92' ^ 5 ' 2 ^ 

which corresponds to forward algorithm 4D 1 £. These are special cases of the general minimal \eyTv\> velocity- type 
algorithm given by by t\ — 0, t$ = 1/{N — 1), 

1 (N-l) 11 , ro , 

v-i = vm = , Vi = — ; r, with evTV = ; r- 5.3 

1 N 2N' N(N-2)' V 2AN{N-2) V ; 

This arbitrary algorithm can serve as a useful check for any general fourth-order, velocity-type algorithm. 
Alternatively, for N — 4, we can allow t 2 to be a free parameter so that 



Theorem 2a then fixes C 2 and evTV with 



and 



U = h, t 3 = 1 - 2t 2 . (5.4) 
c^ = &t 2 {l-t 2 ) 2 (5.5) 



V2 = V3= 12t 2 (Lt 2 y Vi=«4 = 5-«2 (5.6) 

One recognizes that this is the one-parameter algorithm 4BDA first found in Ref. 14 using symbolic algebra. For 
t 2 = 1/2, one recovers integrator 4A; for t 2 = 1/3, one gets back 4D. The advantage of using a variable t 2 is that 
one can use it to minimize the resulting fourth-order error (oftentime to zero) in any specific application. All the 
seven-stage, forward integrators in the velocity form described by Omelyan, Mryglod and Folk (OMF) 1 ^ correspond 
to different ways of choosing t 2 and distributing the commutator term in 4BDA. 

For N — 5, again using t 2 as a parameter, we have t\ = 0, t$ = t 2 , t± = £3 = 1/2 — t 2 , 14.511 with 

V5 = Vi, Va = v 2 , V3 = 1 — 2(vi + v 2 ), and 

v 1 = ^ + C 2 (l-t 2 ), v 2 = -hj 2 . (5.8) 

This is a new one-parameter family of fourth order algorithms with 9 stages or operators. 

To generate position-type algorithms, one can apply Theorem 2b. For N = 3, with v% = 0, v± = v 2 = 1/2, we have 

t 1= i 3 = -(l-— ), *2 = ^=, and e VTV = - — (l--V3). (5.9) 

This produces forward algorithm 4B2LSS corresponding to t 2 = (1 — l/v / 3)/2 in 4BDA. Again, this is a special case 
of the general fourth-order, minimal |e\/Tv| algorithm with v% = 0, Uj = 1/(N — 1), 

(5.10) 




y/N(N-2)- 
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and 

1 / ^N(N-2) \ 
e VTV = -^\l- {N _ 1) j. (5.11) 

For N = 4, Wi = Q and v 2 as the free parameter, invoking Theorem 2b gives 

va = v 2 , V3 = 1 — 2u 2 , (5-12) 

i 2 =i 3 = — 1=, h=t A = \-t 2 (5.13) 
2y/bv 2 2 

and 

evTV = -jo [l-(l-«2)V6t^] ( 5 - 14 ) 

For v 2 = 1/6 and u 2 = 3/8, this reproduces algorithm 4A and 4C— respectively. One again recognizes that the 
above is the one-parameter algorithm 4ACB first derived in Refii^, but now with a much simpler parametrization. 
Algorithm 4ACB covers all the seven-stage, forward fourth-order position-type integrators described by OMF^i. 

For TV = 5, with v 2 as a free parameter, we have v± = 0, v$ = v 2l v 4 — v 3 = 1/2 — v 2l and Theorem 2b produces 
another 9-stage fourth-order algorithm with 



0' = Vl5/!6- 3(v 2 - 1/4) 2 . (5.15) 

h = h, t 4 = t 2 , *3 = 1 - 2(£i + t 2 ), and 

h = \ + C 2 {l-v 2 ), t 2 = ~C 2 . (5.16) 

For N < 5, we have shown above that all fourth-order algorithms are necessarily extended-linear. For N > 5, this 
is not necessary the case. Nevertheless we find that, remarkably, most known N = 5 (9 stages) forward algorithms 
are very close to being extended-linear. For velocity-type, N — 5 extended-linear algorithms, v\ and v 2 are functions 
of t 2 fixed by l|5.8|l . In Fig.l, we compare this predicted relationship with the actual values of vi, v 2 and t 2 of five 
forward, velocity-type, fourth-order algorithms found by OMF— . These are their Eqs.(52)-(56), with their 9, *&, X 
corresponds to t 2 V\, and v 2 respectively. Four out their five algorithms, with v\ in particular, are well described by 
(l5~31) . 

In Fig. 2, we compare the coefficients of all three of OMF's forward, position-type algorithms, Eq.(59)-(61), with 
(|5.16|l . which fixes t\, t 2 as a function of v 2 . Here, their parameters A, p, 8 correspond to v 2 , t\, t 2 respectively. Again, 
t\ is particularly well predicted by l|5.16|l . 

For 11-stage algorithms with N — 6, we have two free parameters t 2 , for velocity type algorithms with 

(j> = 1 - 2*| - 2t\ - (1 - 2t 2 - 2t 3 ) 3 (5.17) 
and two free parameters v 2 , V3 for position type algorithms with 



(j)' = Jl - 2v% - 2w| - (1 - 2v 2 - 2v 3 ) 3 . (5.18) 

Once (j) and (j)' are known, we can determine v\ and v 2 in the case of velocity-type algorithms and t\ and t 2 in the 
case of position-type algorithms. There is one 11-stage velocity algorithm with positive coefficients found by OMF; 
their Eq.(68) with p(= t 2 ) = 0.2029, 9{= t 3 ) = 0.1926, 

i?(= Vi) = 0.0667, and A(= v 2 ) = 0.2620. (5.19) 

The last two values are to be compare with the values given by Theorem 2a below at the same values of t 2 and £3, 

vi = 0.0848, and v 2 = 0.2060. (5.20) 

For OMF's 11-stage, position-type algorithm Eq.(78), with i9(= v 2 ) = 0.1518, A(= v 3 ) = 0.2158, 

p{= h) = 0.0642, and 6{= t 2 ) = 0.1920, (5.21) 

For the same values of v 2 and v 3 , Theorem 2b gives 

t x = 0.0659, and t 2 = 0.1881. (5.22) 

It is remarkable that these 11-stage, fourth-order algorithms derived by complex symbolic algebra, remained very 
close to the values predicted by our extended-linear algorithms. 
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VI. THE STRUCTURE OF NON-FORWARD INTEGRATORS 



Theorem 3 can be used to derive two distinct families of non-forward, fourth-order algorithms. Consider first the 
case of N = 4. For t\ = with symmetric coefficients i 4 = t 2 , the constriants 

2t 2 + h = l (6.1) 

2t\ + 4 = (6.2) 



have unique solutions 



Ea. (|4.14|l then yields 



1 , 2*/3 



and t 3 = -- — — j-. (6.3) 



2 - 2V3 J 2 - 2 1 / 3 



1 1 l(2 1 /3_i) 
Vl=Vi= 22—¥T^ v * = v * = — 2 2 _ 2 i/3 ■ ( 6 - 4 ) 

One recognizes that we have just derived the well known fourth-order Forest-Ruth integrator—. Note that there is 
complete symmetry between {U} and {vi}. For position type algorithm, we simply interchange the values of and 

Vi. 

There are no symmetric solutions for N = 5, for the same reason that there are also no solutions for N — 3. For 
N = 2k, we have the general condition 



2 tj +**+! = !, (6.5) 



i=2 



2^^+4+1 = 0, (6.6) 



i=2 

which can be solved by introducing real parameters on for i = 2 to k with a2 = 1, 

U = a { t 2 , (6.7) 

so that 

1/3 



* fc+ i = -2 1/3 (j>?) t a , 



(6.8) 



h = jyj. (6.9) 

2 (El 2 - 2V3 (J2l 2 a ? 



These solutions generalize the fourth-order Forest-Ruth integrator to arbitrary N. 
For N = 2k + 1, k > 2, again introducing 1)6.7(1 for i = 2 to k with = 1, we have 



\ " .3 



i fe+ i = - *2, (6.10) 



\i=2 



, a 3 



V3 - 



2 (EL2 a *) - 2 (E 8 fe = 2 "i 

This is a new class of fourth-order algorithm possible only for N odd and greater than five. 



(6.11) 
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VII. CONCLUSIONS 

Most of the machinery for tracking coefficients were developed in Rcf. 18 for the purpose of providing a constructive 
proof of the Sheng-Suzuki theorem. The advantage of this constructive approach is that we can obtain explicit lower 
bounds on the the second-order error coefficients. Here, by imposing the extended-linear relationship between {ti} 
and {vi}, these bounds become the actual error coefficients and provide a complete characterization all fourth-order 
symplectic integrators for arbitrary number of operators. The most satisfying aspect of this work is that most fourth 
order integrators can now be derived analytically without recourse to symbolic algebra or numerical root-finding. We 
have also provided explicit construction of many new classes of fourth-order algorithms. 

For N = 5, 6, corresponding to 9 and 11 operators, we have shown that many fourth-order algorithms found by 
Omelyan, Mryglod and Folkii are surprisely close to the predicted coefficients of our theory, suggesting that the 
extended-linear relation between coefficients may be the dominate solution of the order-condition. 

The expansion (|3.2|l may hold similar promise for characterizing sixth order algorithms by introducing extended- 
quadratic or higher order relationships between the two sets of coefficients. 
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FIG. 1: Comparing the coefficients of five, 9-stage, velocity-type, fourth-order forward integrators of Omelyan, Mryglod and 
Folk— (filled circles and squares) , with the analytical prediction of extended-linear symplectic integrators (solid lines) . 
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FIG. 2: Comparing the coefficients of three, 9-stage, position-type, fourth-order forward integrators of Omelyan, Mryglod and 
Folk— (filled circles and squares) , with the analytical prediction of extended-linear symplectic integrators (solid lines) . 



